Jean-Christophe Loiseau
Maître de Conférences
📧 : jean-christophe.loiseau@ensam.eu
📅 : December 19^th 2019

# --> Julia packages for the plots.
using Plots
using LaTeXStrings

# --> Julia package to integrate forward in time ODEs.
using DifferentialEquations

# -->
using DynamicalSystems

Introduction to nonlinear physics, dynamical systems and chaos theory

The aim of this notebook, in association with the first lecture of the dynamical system class I teach at Arts et Métiers, is to illustrate some of the various concepts that will be taught during the course on a set of relatively simple, yet illustrative, examples. For that purpose, the following systems will be considered:

• Double pendulum: canonical example from classical physics belonging to the class of Hamiltonian systems, i.e. systems for which a quantity (namely the Hamiltonian, here the sum of the potential and kinetic energies) is conserved.
• Lorenz system: canonical example from dynamical system theory belonging to the class of dissipative systems of particular historical importance.
• Logistic map: Another classical example from dynamical system theory used to illustrate the concept of bifurcations. Note that, contrary to all other examples considered herein, this one is formulated in a discrete-time framework.

Provided you have installed julia and the required packages (i.e. Plots.jl, PyPlot.jl, DifferentialEquations.jl and DynamicalSystems.jl), this notebook should be self-contained. All of the pieces of code needed to run the various examples are included. If you encounter any problem when running this notebook, please do not hesitate to get back to me either by email (📧) or directly on GitHub and I'll try to fix it.

License: This work is distributed under the Creative Commons BY-NC-SA 4.0 license.

---

The double pendulum

The double pendulum is a canonical example from classical mechanics which can be used to model various mechanical systems. It belongs to the class of so-called Hamiltonian systems for which a quantity (namely the Hamiltonian, here the sum of kinetic and potential energies) is conserved. Despite its apparent simplicity, this system can nonetheless exhibit extremely complex dynamics. Defining the angular positions $\theta_1$ and $\theta_2$ as shown in the left figure and the corresponding angular velocities $\dot{\theta}_1$ and $\dot{\theta}_2$, the Lagrangian $\mathcal{L}(\theta_1, \theta_2, \dot{\theta}_1, \dot{\theta}_2)$ of this system is the difference of the kinetic energy given by

$$ \mathcal{T} = \frac{1}{2} m_1 l_1^2 \dot{\theta}_1^2 + \frac{1}{2} m_2 \left( l_1^2 \dot{\theta}_1^2 + l_2^2 \dot{\theta}_2^2 + 2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \cos \left( \theta_1 - \theta_2 \right) \right) $$

and potential energy

$$ \mathcal{V} = (m_1 + m_2) g l_1 \cos \theta_1 + m_2 g l_2 \cos \theta_2 $$

such that

$$ \mathcal{L} = \mathcal{T} - \mathcal{V}. $$

The equations of motion are then given by

$$ \frac{\mathrm{d}}{\mathrm{d}t} \left( \frac{\partial \mathcal{L}}{\partial \dot{\theta}_i} \right) - \frac{\partial \mathcal{L}}{\partial \theta_i} = 0. $$

Similarly, the equations of motions could be derived from the Hamiltonian $\mathcal{H}(p_1, p_2, \theta_1, \theta_2)$

$$ \mathcal{H} = \frac{m_2 l_2^2 p_1^2 + (m_1 + m_2) p_2^2 - 2 m_2 l_1 l_2 p_1 p_2 \cos(\theta_1 - \theta_2)}{2 l_1^2 l_2^2 m_2 (m_1 + m_2 \sin^2(\theta_1 - \theta_2)} - (m_1 + m_2) g l_1 \cos \theta_1 - m_2 g l_2 \cos \theta_2 $$

with

$$ p_i = \frac{\partial \mathcal{L}}{\partial \dot{\theta}_i} $$

the generalized momenta. The equations of motion are then given by

$$ \begin{aligned} & \frac{\mathrm{d} \theta_i}{\mathrm{d}t} = \frac{\partial \mathcal{H}}{\partial p_i} \\ & \frac{\mathrm{d} p_i}{\mathrm{d}t} = - \frac{\partial \mathcal{H}}{\partial \theta_i}. \end{aligned} $$

This Hamiltonian formulation of the problem will be used in julia below to simulate this system. The next few cells set up the problem and simulate the evolution of the system for various initial conditions.

# --> Define the mass of each pendulum.
m₁ = 1.0
m₂ = 1.0

# --> Define the length of each pendulum.
l₁ = 1.0
l₂ = 1.0

# --> Intensity of the gravity.
g = 9.81

# --> Pack all the parameters in a tuple for the sake of simplicity.
params = (m₁, m₂, l₁, l₂, g)

(1.0, 1.0, 1.0, 1.0, 9.81)

function Hamiltonian_double_pendulum(p, q, params)
    
    """
    This function computes the Hamiltonian for the double pendulum.
    
    INPUTS
    ------
    
    p :
    
    q : Angular positions of the two pendulums.
    
    params : Parameters defining the system.
    
    OUTPUT
    ------
    
    H : Scalar
        Hamiltonian of the system.
    """
   
    # --> Unpack parameters.
    m₁, m₂, l₁, l₂, g = params
    
    # --> Compute the Hamiltonian.
    H₁ = m₂ * l₂^2 * p[1].^2 + (m₁ + m₂) * l₁^2 * p[2].^2 - 2m₂ * l₁ * l₂ * p[1] * p[2] * cos(q[1] - q[2])
    H₂ = 2m₂ * l₁^2 * l₂^2 * (m₁ + m₂ * sin(q[1] - q[2])^2)
    H₃ = -(m₁ + m₂) * g * l₁ * cos(q[1]) - m₂ * g * l₂ * cos(q[2])
    
    H = H₁ / H₂ + H₃
    
    return H
end

Hamiltonian_double_pendulum (generic function with 1 method)

# --> Temporal window over which the system is simulated.
tspan = (0.0, 25.0)

# --> Instants of time at which the evolution is sampled.
Δt = 0.025;
t = collect(tspan[1]:Δt:tspan[2]);

# --> Zero initial condition for the angular velocities.
#     NOTE: Same condition for all simulations.
p₀ = [0.0; 0.0]

# --> Small amplitudes simulation.
q₀ = [π/5 ; π/5]

prob₀ = HamiltonianProblem(
    Hamiltonian_double_pendulum, 
    p₀, q₀, 
    tspan, 
    p=params
)

sol₀ = solve(prob₀, McAte5(), tol=1e-8, dt=1e-3);

# --> Intermediate amplitudes simulation.
q₁ = [2π/5 ; 2π/5]

prob₁ = HamiltonianProblem(
    Hamiltonian_double_pendulum, 
    p₀, q₁, 
    tspan, 
    p=params
)

sol₁ = solve(prob₁, McAte5(), tol=1e-8, dt=1e-3);

# --> Intermediate amplitudes simulation (bis).
q₂ = [3π/5 ; 3π/5]

prob₂ = HamiltonianProblem(
    Hamiltonian_double_pendulum, 
    p₀, q₂, 
    tspan, 
    p=params
)

sol₂ = solve(prob₂, McAte5(), tol=1e-8, dt=1e-3);

# --> Large amplitudes simulation
q₃ = [4π/5 ; 4π/5]

prob₃ = HamiltonianProblem(
    Hamiltonian_double_pendulum, 
    p₀, q₃, 
    tspan, 
    p=params
)

sol₃ = solve(prob₃, McAte5(), tol=1e-8, dt=1e-3);

function polar_to_cartesian(t, sol, params)
    
    """
    Utility function to go from polar to cartesian coordinates.
    """
    
    # --> Unpack parameters.
    m₁, m₂, l₁, l₂, g = params

    # -->
    x₁ = l₁ * sin.(sol(t)[3, :])
    y₁ = -l₁ * cos.(sol(t)[3, :])
    
    # -->
    x₂ = x₁ + l₂ * sin.(sol(t)[4, :])
    y₂ = y₁ - l₂ * cos.(sol(t)[4, :])
    
    return [x₁ y₁ x₂ y₂]
end

polar_to_cartesian (generic function with 1 method)

# --> Transform each simulation into the cartesian reference frame.
x₀ = polar_to_cartesian(t, sol₀, params)
x₁ = polar_to_cartesian(t, sol₁, params)
x₂ = polar_to_cartesian(t, sol₂, params)
x₃ = polar_to_cartesian(t, sol₃, params);

function plot_pendulum(i, x)

    """
    Utility function to plot the evolution of the double pendulum.
    """
    
    # --> Plot the position of the first pendulum.
    p = plot(
        [0, x[i, 1]], [0, x[i, 2]],
        size=(512, 512),
        xlim=(-2.2, 2.2),
        ylim=(-2.2, 2.2),
        markersize=15,
        markershape=:circle,
        label="",
        framestyle=:none,
    )
    
    # --> Plot the position of the second pendulum.
    p = plot!(
        [x[i, 1], x[i, 3]], [x[i, 2], x[i, 4]],
        markersize=15,
        markershape=:circle,
        label="",
        aspect_ratio=:equal,
    )
    
    # --> Plot the streak line effect.
    if i > 9*2
        p = plot!([x[i-3*2:i, 3]], [x[i-3*2:i, 4]],alpha = 0.15,linewidth = 2, color = :red,label ="");
        p = plot!([x[i-5*2:i-3*2, 3]], [x[i-5*2:i-3*2, 4]],alpha = 0.08,linewidth = 2, color = :red,label ="");
        p = plot!([x[i-7*2:i-5*2, 3]], [x[i-7*2:i-5*2, 4]],alpha = 0.04,linewidth = 2, color = :red, label ="");
        p = plot!([x[i-9*2:i-7*2, 3]], [x[i-9*2:i-7*2, 4]],alpha = 0.01,linewidth = 2, color = :red, label="");
    end
    
    return p
end

plot_pendulum (generic function with 1 method)

gr()

# --> Loop through the frames.
anim = @animate for i = 1:length(t)
    # --> Plot each simulation.
    p₀ = plot_pendulum(i, x₀);
    p₁ = plot_pendulum(i, x₁);
    p₂ = plot_pendulum(i, x₂);
    p₃ = plot_pendulum(i, x₃);

    p = plot(
        p₀, p₁, p₂, p₃, 
        layout=(2, 2), 
        size=(1024, 1024)
    );
end

# --> Display the gif animation.
gif(anim, fps=30)

┌ Info: Saved animation to 
│   fn = /home/loiseau/Documents/Teaching/Nonlinear_physics/Lecture_1/Scripts/tmp.gif
└ @ Plots /home/loiseau/.julia/packages/Plots/qZHsp/src/animation.jl:98

The animation above shows the evolution of the exact same double pendulum (i.e. same masses and same lengths) for four different initial conditions. While the upper two evolutions appear to be quite periodic, the lower two show a more erratic behaviour hardly predictable despite the deterministic nature of the equations of motion. As we will see later on during the course, the lower two pendula exhibit chaotic dynamics.

Sensitivity to initial conditions

One of the hallmarks of chaos is the so-called sensitivity to initial conditions: if a system is chaotic, even the smallest difference in the initial conditions can lead to vastly different evolutions. This is illustrated below where two double pendula are released at angles differing by less than $\frac{\pi}{1000}$.

# --> Temporal window over which the system is simulated.
tspan = (0.0, 25.0)

# --> Perturbation.
ϵ = 1e-3

# --> Instants of time at which the evolution is sampled.
Δt = 0.025;
t = collect(tspan[1]:Δt:tspan[2]);

# --> Zero initial condition for the angular velocities.
#     NOTE: Same condition for all simulations.
p₀ = [0.0 ; 0.0]

# --> Reference solution.
q₀ = [3π/5 ; 3π/5]

prob₀ = HamiltonianProblem(
    Hamiltonian_double_pendulum, 
    p₀, q₀, 
    tspan, 
    p=params
)

sol₀ = solve(prob₀, McAte5(), tol=1e-8, dt=1e-3);

# --> Perturbed initial condition.
q₁ = q₀ .+ ϵ

prob₁ = HamiltonianProblem(
    Hamiltonian_double_pendulum, 
    p₀, q₁, 
    tspan, 
    p=params
)

sol₁ = solve(prob₁, McAte5(), tol=1e-8, dt=1e-3);

# --> Transform each simulation into the cartesian reference frame.
x₀ = polar_to_cartesian(t, sol₀, params);
x₁ = polar_to_cartesian(t, sol₁, params);

gr()

# --> Loop through the frames.
anim = @animate for i = 1:length(t)
    # --> Plot each simulation.
    p₀ = plot_pendulum(i, x₀);
    p₁ = plot_pendulum(i, x₁);

    p = plot(
        p₀, p₁, 
        layout=(1, 2), 
        size=(1024, 512)
    );
end

# --> Display the gif animation.
gif(anim, fps=30)

┌ Info: Saved animation to 
│   fn = /home/loiseau/Documents/Teaching/Nonlinear_physics/Lecture_1/Scripts/tmp.gif
└ @ Plots /home/loiseau/.julia/packages/Plots/qZHsp/src/animation.jl:98

As can be observed, although the two pendula initially follow the same trajectory, their evolutions quickly differ after only a few oscillations. This is a striking illustration of sensitivity to initial conditions exhibited by all chaotic dynamical systems.

---

Lorenz system

The second dynamical system we'll look at is the Lorenz system. As we will show in a coming lesson, this system can be derived from the Navier-Stokes equations governing the motion of the fluid between two parallel plates with the lower one being heated up while the upper one is cooled down. The governing equations read

$$ \begin{aligned} & \frac{\mathrm{d}x}{\mathrm{d}t} = \sigma (y - x) \\ & \frac{\mathrm{d}y}{\mathrm{d}t} = x (\rho - z) - y \\ & \frac{\mathrm{d}z}{\mathrm{d}t} = xy - \beta z \end{aligned} $$

where $x$, $y$ and $z$ are the three state variables and $\sigma$, $\rho$ and $\beta$ are user-defined parameters. In what follows, we will consider the canonical set of parameters given by $\sigma = 10$, $\rho = 28$ and $\beta = \frac{8}{3}$. It should be noted that, contrary to the double pendulum system, the Lorenz system is not a Hamiltonian system but a dissipative system, i.e. energy is not conserved. The next few cells set up the problem and simulate it for random initial conditions.

function lorenz!(du, u, p, t)
   
    # --> Unpack parameters.
    σ, ρ, β = p
    
    # --> Governing equations.
    du[1] = σ * (u[2] - u[1])
    du[2] = u[1] * (ρ - u[3]) - u[2]
    du[3] = u[1] * u[2] - β * u[3]
    
    return du
end

lorenz! (generic function with 1 method)

# --> Setup the parameters.
σ, ρ, β = 10.0, 28.0, 8.0/3.0
params = (σ, ρ, β)
tspan = (0.0, 50.0)

# --> Initial condition.
u₀ = [10.0 ; 10.0 ; 10.0]

# --> Setup and solve the problem.
prob = ODEProblem(lorenz!, u₀, tspan, params)
sol = solve(prob);

Let us now plot the evolution of each state variable.

gr()

# --> Plot the time-series of the three state variables.
p₀ = plot(
    sol, vars=(0, 1), 
    legend=:none,
    ylabel=L"x(t)"
)
p₁ = plot(
    sol, vars=(0, 2),
    ylabel=L"y(t)",
    legend=:none
)
p₂ = plot(
    sol, vars=(0, 3), 
    ylabel=L"z(t)",
    legend=:none
)

p = plot(
    p₀, p₁, p₂,
    layout=(3, 1),
    xlim=(0, 50),
    size=(800, 400)
)

Although such a visualization is easy to do, a better understanding of the system's dynamics can be provided by plotting its dynamics in the phase space as shown below.

gr()

p₀ = plot(
    sol, vars=(1, 2),
    legend=:none,
    xlabel=L"x",
    ylabel=L"y",
)

p₁ = plot(
    sol, vars=(1, 3),
    legend=:none,
    xlabel=L"x",
    ylabel=L"z",
)

p₂ = plot(
    sol, vars=(2, 3),
    legend=:none,
    xlabel=L"y",
    ylabel=L"z",
)

p = plot(
    p₀, p₁, p₂,
    layout=(1, 3),
    size=(900, 300)
)

Quite clearly, the oscillations observed in the time-series of the various state variables correspond to oscillations of the system between two different states. As we will investigate during the dedicated lesson, the apparently random switches between the two wings of this so-called strange attractor are not actually random.

Finally, let us consider an ensemble of random initial conditions and see their evolution. This is done below by setting up a Monte-Carlo problem in julia.

function prob_func(prob, i, repeat)
    """Utility function to generate random initial condition"""
    @. prob.u0 = 40 * [rand()-0.5, rand()-0.5, rand()]
    return prob
end

# --> Define and solve the problem for an ensemble of random initial conditions.
ensemble_prob = EnsembleProblem(prob, prob_func=prob_func)
sim = solve(ensemble_prob, Tsit5(), trajectories=50);

gr()

# --> Loop through the frames
Δt = 0.025
n = Int(tspan[2] / Δt)

anim = @animate for i = 100:n
    
    p = plot(
        sim, vars=(1, 3), tspan=((i-20)*Δt, i*Δt),
        xlabel=L"x",
        ylabel=L"z",
        xlim=(-20, 20),
        ylim=(0, 50),
        legend=:none,
        size=(384, 384),
        seriescolor=:auto,
        color_palette=:blues
    )
    
end

# --> Display the gif animation.
gif(anim, fps=30)

┌ Info: Saved animation to 
│   fn = /home/loiseau/Documents/Teaching/Nonlinear_physics/Lecture_1/Scripts/tmp.gif
└ @ Plots /home/loiseau/.julia/packages/Plots/qZHsp/src/animation.jl:98

After some initial transients, all of the simulations converge to the same attractor. This attractor has peculiar properties. The most suprising one is that its dimension is neither 2 (not a surface) nor 3 (not a volume) but in-between. To be precise, its dimension is 2.06. Non-integer dimensions are another hallmarks of chaotic dynamics closely related to the concept of fractal geometry.

---

Logistic map

The last system we will look at is the logistic map, another historically relevant dynamical system closely related to the Mandelbrot set. Contrary to all of the other systems considered herein that are governed by equations of the form $\dot{x} = f(x)$, the logitic map does not evolve continuously in time. Instead, it evolves discretely. It is a discrete-time nonlinear system whose evolution is governed by

$$ x_{k+1} = \mu x_k ( 1 - x_k) $$

where $\mu \in \left[0, 4 \right[$ is the carying capacity and $x_k \in \left[0, 1 \right]$ is the population at time $t_k$. As we will see throughout the course, the logistic equation has played a crucial role in the establishement of the chaos theory precisely due to its simplicity. The next few cells set up the problem and simulate it for various values of the carying capacity $\mu$.

# -->
logistic_map(x, p, t) = p[1] * x * (1.0 - x)

# -->
state, p = 0.01, [0.1]
dynsys = DiscreteDynamicalSystem(logistic_map, state, p)

1-dimensional discrete dynamical system
 state:       0.01
 e.o.m.:      logistic_map
 in-place?    false
 jacobian:    ForwardDiff
 parameters:  [0.1]

function plot_bifurcation_diagram(dynsys, μ_range)
    
    # -->
    n, Ttr = 2000, 20000
    
    # -->
    output = orbitdiagram(dynsys, 1, 1, μ_range, n=n, Ttr=Ttr)
    
    # -->
    l = length(μ_range)
    μ = Vector{Float64}(undef, n*l)
    x = copy(μ)
    for j = 1:l
        μ[(1 + (j-1)*n):j*n] .= μ_range[j]
        x[(1 + (j-1)*n):j*n] .= output[j]
    end
    
    pyplot()
    p = scatter(
        μ, x,
        xlim=(minimum(μ), maximum(μ)), xlabel=L"\mu", tickfontsize=12,
        ylim=(0, 1), ylabel=L"\bar{x}",
        markershape=:circle , markersize=0.002, markercolor=1, markerstrokecolor=1, alpha=0.001,
        legend=:none, framestyle=:box, size=(800, 300)
    )
    
    return p
end

plot_bifurcation_diagram (generic function with 1 method)

plot_bifurcation_diagram(dynsys, collect(1:0.001:3.999))

This plot is known as the bifurcation diagram or orbit diagram of the logistic map. The horizontal axis corresponds to the different values of $\mu$ that have been considered while the asymptotic state of the system (after some transients have faded) are reported on the vertical axis. For $1 \leq \mu \leq 3$, it can be seen that the system is characterized by a steady state (i.e. only a single branch in this bifurcation diagram). For $3 \leq \mu \leq 1 + \sqrt{6}$, the system switches periodically between two different states. For this range of $\mu$, the system exhibits some period-2 dynamics. This change of dynamics, between steady state for $\mu \leq 3$ to period-2 dynamics for $3 \leq \mu$ is known as a bifurcation and $\mu = 3$ is called a bifurcation point. As $\mu$ increases, the system experiences a whole sequence of other bifurcations, yielding it to exhibit increasingly complex dynamics despite the simplicity of the model.

It must be noted finally that this bifurcation diagram exhibits some kind of self-similarity or fractal structure as illustrated by the two magnifications of the bifurcation diagram below. Each branch of this diagram can be understood as a miniaturized version of the whole diagram itself. This property, common to a number of different dynamical systems, can be characterized by the so-called Feigenbaum constants which are just as universal as $\pi$ or $e$. Computing these constants numerically will be the subject of one of the numerical projects proposed during the class.

p = plot_bifurcation_diagram(dynsys, collect(3:0.0001:3.6))
p = plot!(ylim=(2/3, 0.9))

p = plot_bifurcation_diagram(dynsys, collect(1 + √6:0.00005:3.6))
p = plot!(ylim=(0.85, 0.9))